# load required libraries
library("lme4")
library("ggplot2")
library("here")
library("tidyverse")
library("emmeans")
library("brms")
library("DHARMa")
library("glmmTMB")
library("performance")
library("scales")

1 Exploration and visualisation

A pseudonymized/codified dataset was used as a precaution due to ethical considerations (patient-level clinical data). The full dataset is available upon request (see manuscript data availability statement).

# load the data
df <- read_csv(here("./data/Ig_detection_codified.csv"), show_col_types = FALSE) %>%
  # rename(Experiment = `Experiment (= Sample collection Number)`) %>% 
  rename(HBB.genotype = `HBB Genotype`,
         `Raw IgM+ RBC counts_IgM plate` = `Raw IgG+ RBC counts_IgM plate`) %>% 
  mutate(across(c(`Date Flow cytometry experiment`), as.Date)) %>% 
  select(
    ID, HBB.genotype, Sample_type,
    `Raw IgG+ RBC counts_IgG plate`,    `Raw IgG+ iRBC counts`, `Raw IgG+ Asexual iRBC counts`, `Raw IgG+ Sexual iRBC counts`,
    `Raw IgM+ RBC counts_IgM plate`,    `Raw IgM+ iRBC counts`, `Raw IgM+ Asexual iRBC counts`, `Raw IgM+ Sexual iRBC counts`,
    `Raw IgG- iRBC counts`, `Raw IgG- Asexual iRBC counts`, `Raw IgG- Sexual iRBC counts`,  
    `Raw IgM- iRBC counts`, `Raw IgM- Asexual iRBC counts`, `Raw IgM- Sexual iRBC counts`,
    Age_group, Gender, Village, CollectionSite, 
    Pf_infection, Parasites_µL, Gametocytes_µL, iRBCs_IgG_PercentageCells, iRBCs_IgM_PercentageCells
  ) %>% 
  rename(
    raw_igGpos_RBC = `Raw IgG+ RBC counts_IgG plate`,
    raw_igGpos_iRBC = `Raw IgG+ iRBC counts`,
    raw_igGpos_asexual = `Raw IgG+ Asexual iRBC counts`,
    raw_igGpos_sexual = `Raw IgG+ Sexual iRBC counts`,
    raw_igMpos_RBC = `Raw IgM+ RBC counts_IgM plate`,
    raw_igMpos_iRBC = `Raw IgM+ iRBC counts`,
    raw_igMpos_asexual = `Raw IgM+ Asexual iRBC counts`,
    raw_igMpos_sexual = `Raw IgM+ Sexual iRBC counts`,
    raw_igGneg_iRBC = `Raw IgG- iRBC counts`,   
    raw_igGneg_asexual = `Raw IgG- Asexual iRBC counts`,
    raw_igGneg_sexual = `Raw IgG- Sexual iRBC counts`,
    raw_igMneg_iRBC = `Raw IgM- iRBC counts`,
    raw_igMneg_asexual = `Raw IgM- Asexual iRBC counts`,
    raw_igMneg_sexual = `Raw IgM- Sexual iRBC counts`,
  ) %>% 
  mutate(HBB.genotype = case_when(
    HBB.genotype == 0 ~ "HbAA",
    HBB.genotype == 1 ~ "HbAC",
    HBB.genotype == 2 ~ "HbAS",
    HBB.genotype == 3 ~ "HbCC",
    HBB.genotype == 4 ~ "HbSC",
  )) %>%
  mutate(Pf_infection = case_when(
    Pf_infection == 0 ~ "NO",
    Pf_infection == 1 ~ "YES",
  )) %>%
  mutate(across(c(ID, HBB.genotype, Sample_type, Gender, Village,
                  CollectionSite, Pf_infection, Age_group), 
                as.factor)) %>%
  mutate(HBB.genotype = fct_relevel(HBB.genotype, "HbAA")) %>% # set HbAA as baseline reference
  drop_na(starts_with("raw_"))
  
all(df$HBB.genotype == df$Sample_type)

1.1 IgG

plot_func <-
  function(data, x, y, title){
    data %>%
      ggplot(
        aes(x = {{x}},
            y = {{y}}))+
      geom_boxplot(aes(fill = {{x}}), outlier.size = 2, outlier.colour = "grey",
                   width = 0.8, position = position_dodge(width = 0.8)) +  # Reduce space between boxplots
      geom_point(aes(fill = {{x}}), alpha = 0.8, shape = 21, colour= "#115e67a6",
                 position = position_jitter(width = 0.1, seed = 42)
                 ) +
      scale_fill_manual(values = 
                          c("HbAA" ="paleturquoise","AA" ="paleturquoise", 
                            "HbAC" = "palegreen", 
                            "HbAS" = "thistle",
                            "HbCC" = "lightpink1",
                            "HbSC"= "peachpuff")) +
      scale_y_continuous(
        # limits = c(0, 1), 
        labels = label_percent()) +
      theme_minimal() +
      theme(
        axis.text.x = element_text(angle = 45, hjust = 1, size = 14),  # Increase font size for x-axis ticks
        axis.text.y = element_text(size = 14),  # Increase font size for y-axis ticks
        axis.title.x = element_text(size = 16),  # Increase font size for x-axis title
        axis.title.y = element_text(size = 16),  # Increase font size for y-axis title
        strip.text = element_text(size = 12),
        panel.grid.major = element_blank(),   # Remove major gridlines
        panel.grid.minor = element_blank(),   # Remove minor gridlines
        axis.line = element_line(color = "black")  # Add black axis lines
      ) +
      labs(title = title)
}

plot_func(data = df, x = HBB.genotype, y = raw_igGpos_iRBC / (raw_igGpos_iRBC + raw_igGneg_iRBC), "iRBC IgG+ / iRBCs")

plot_func(data = df, x = HBB.genotype, y = raw_igGpos_asexual / (raw_igGpos_asexual + raw_igGneg_asexual), "Asexual IgG+ / asexual iRBCs")

plot_func(data = df, x = HBB.genotype, y = raw_igGpos_sexual / (raw_igGpos_sexual + raw_igGneg_sexual), "Sexual iRBC IgG+ / sexual iRBCs")

1.2 IgM

plot_func(data = df, x = HBB.genotype, y = raw_igMpos_iRBC / (raw_igMpos_iRBC + raw_igMneg_iRBC), "iRBC IgM+ / iRBCs")

plot_func(data = df, x = HBB.genotype, y = raw_igMpos_asexual / (raw_igMpos_asexual + raw_igMneg_asexual), "Asexual IgM+ / asexual iRBCs")

plot_func(data = df, x = HBB.genotype, y = raw_igMpos_sexual / (raw_igMpos_sexual + raw_igMneg_sexual), "Sexual iRBC IgM+ / sexual iRBCs")

1.3 Pf infection

ggplot(df,
       aes(x = Pf_infection, y = raw_igGpos_iRBC / (raw_igGpos_iRBC + raw_igGneg_iRBC))) +
      geom_boxplot(width = 0.8, position = position_dodge(width = 0.8)) + 
      geom_point(alpha = 0.8, shape = 21, colour= "#115e67a6",
                 position = position_jitter(width = 0.1, seed = 42))

ggplot(df,
      aes(x = Pf_infection, y = iRBCs_IgG_PercentageCells/100)) +
      geom_boxplot(width = 0.8, position = position_dodge(width = 0.8)) + 
      geom_point(alpha = 0.8, shape = 21, colour= "#115e67a6",
                 position = position_jitter(width = 0.1, seed = 42))

ggplot(df,
       aes(x = Pf_infection, y = raw_igMpos_iRBC / (raw_igMpos_iRBC + raw_igMneg_iRBC))) +
      geom_boxplot(width = 0.8, position = position_dodge(width = 0.8)) + 
      geom_point(alpha = 0.8, shape = 21, colour= "#115e67a6",
                 position = position_jitter(width = 0.1, seed = 42)) +
        scale_y_continuous(
        limits = c(0, .15),
        labels = label_percent())

ggplot(df,
      aes(x = Pf_infection, y = iRBCs_IgM_PercentageCells/100)) +
      geom_boxplot(width = 0.8, position = position_dodge(width = 0.8)) + 
      geom_point(alpha = 0.8, shape = 21, colour= "#115e67a6",
                 position = position_jitter(width = 0.1, seed = 42)) +
          scale_y_continuous(
        limits = c(0, .15),
        labels = label_percent())

Comparison between percentage and raw values.

plot_func(data = df, x = Pf_infection, y = raw_igGpos_iRBC / (raw_igGpos_iRBC + raw_igGneg_iRBC), "iRBC IgG+ / iRBCs")
## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's fill values.
## No shared levels found between `names(values)` of the manual scale and the
## data's fill values.
## No shared levels found between `names(values)` of the manual scale and the
## data's fill values.

plot_func(data = df, x = Pf_infection, y = iRBCs_IgG_PercentageCells, "iRBCs_IgG_PercentageCells")
## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's fill values.
## No shared levels found between `names(values)` of the manual scale and the
## data's fill values.
## No shared levels found between `names(values)` of the manual scale and the
## data's fill values.

plot_func(data = df, x = Pf_infection, y = raw_igMpos_iRBC / (raw_igMpos_iRBC + raw_igMneg_iRBC), "iRBC IgM+ / iRBCs")
## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's fill values.
## No shared levels found between `names(values)` of the manual scale and the
## data's fill values.
## No shared levels found between `names(values)` of the manual scale and the
## data's fill values.

plot_func(data = df, x = Pf_infection, y = iRBCs_IgM_PercentageCells, "iRBCs_IgM_PercentageCells")
## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's fill values.
## No shared levels found between `names(values)` of the manual scale and the
## data's fill values.
## No shared levels found between `names(values)` of the manual scale and the
## data's fill values.

Figure for publication:

# Color palette
fill_colors <- c("NO"="#FFDAB5",
                 "YES" ="#7BAFD4")
outline_colors <- c(
  "NO"="#FF8C00",
  "YES" ="#375D81")


# Updated outline colors (darker tones)
outline_colors <- c(
  "NO"="#FF8C00",
  "YES"= "#375D81"
)

summary_stats <- df %>%
  group_by(Pf_infection) %>%
  summarise(
    ymin = min(raw_igGpos_iRBC / (raw_igGpos_iRBC + raw_igGneg_iRBC), na.rm = TRUE),
    ymax = max(raw_igGpos_iRBC / (raw_igGpos_iRBC + raw_igGneg_iRBC), na.rm = TRUE),
    ymed = median(raw_igGpos_iRBC / (raw_igGpos_iRBC + raw_igGneg_iRBC), na.rm = TRUE)
  )

p <- ggplot(df, aes(x = Pf_infection , y = raw_igGpos_iRBC / (raw_igGpos_iRBC + raw_igGneg_iRBC) )) +
  
  geom_jitter(aes(color = Pf_infection, fill = Pf_infection),
              width = 0.15, size = 3.5, shape = 21, alpha = 0.9, stroke = 0.8) +
  
  geom_errorbar(data = summary_stats, aes(ymin = ymin, ymax = ymax, x = Pf_infection),
                width = 0.2, color = "grey60", inherit.aes = FALSE) +
  
  geom_point(data = summary_stats, aes(y = ymed, x = Pf_infection),
             shape = 95, size = 8, color = "black", inherit.aes = FALSE) +
  
  scale_fill_manual(values = fill_colors) +
  scale_color_manual(values = outline_colors) + 
  
  # scale_y_continuous(limits = c(0, 30), breaks = seq(0, 30, 10), expand = expansion(mult = c(0.03, 0.03))) +
  scale_y_continuous(
    # limits = c(0, .15),
    labels = label_percent()
    ) +
  
  geom_hline(yintercept = 0, linetype = "dotted", color = "grey40") +
  
  theme_minimal(base_size = 15) +
  theme(
    panel.grid = element_blank(),
    axis.line = element_line(color = "black"),
    axis.text.y = element_text(color = "black", size = 14),  # Bigger Y-axis numbers
    axis.text.x = element_text(color = "black", size = 14),  # Bigger X-axis category names
    axis.title = element_text(color = "black", size = 16),   # Bigger axis titles
    legend.position = "none"
  ) +
  
  labs(x = "Pf infection status", y = "Total IgG+ Troph-iRBCs (% cells)")

p

ggsave(here("./figures/pf_infection_IgG.png"))
## Saving 7 x 5 in image
# Summary stats
summary_stats <- df %>%
  group_by(Pf_infection) %>%
  summarise(
    ymin = min(raw_igMpos_iRBC / (raw_igMpos_iRBC + raw_igMneg_iRBC), na.rm = TRUE),
    ymax = max(raw_igMpos_iRBC / (raw_igMpos_iRBC + raw_igMneg_iRBC), na.rm = TRUE),
    ymed = median(raw_igMpos_iRBC / (raw_igMpos_iRBC + raw_igMneg_iRBC), na.rm = TRUE)
  )

# Color palette
fill_colors <- c("NO"="#FFDAB5",
                 "YES" ="#7BAFD4")
outline_colors <- c(
  "NO"="#FF8C00",
  "YES" ="#375D81")


# Updated outline colors (darker tones)
outline_colors <- c(
  "NO"="#FF8C00",
  "YES"= "#375D81"
)


p <- ggplot(df, aes(x = Pf_infection , y = raw_igMpos_iRBC / (raw_igMpos_iRBC + raw_igMneg_iRBC) )) +
  
  geom_jitter(aes(color = Pf_infection, fill = Pf_infection),
              width = 0.15, size = 3.5, shape = 21, alpha = 0.9, stroke = 0.8) +
  
  geom_errorbar(data = summary_stats, aes(ymin = ymin, ymax = ymax, x = Pf_infection),
                width = 0.2, color = "grey60", inherit.aes = FALSE) +
  
  geom_point(data = summary_stats, aes(y = ymed, x = Pf_infection),
             shape = 95, size = 8, color = "black", inherit.aes = FALSE) +
  
  scale_fill_manual(values = fill_colors) +
  scale_color_manual(values = outline_colors) +
  
  # scale_y_continuous(expand = expansion(mult = c(0.03, 0.03))) +
  scale_y_continuous(
    # limits = c(0, .15),
    labels = label_percent()
    ) +
  
  geom_hline(yintercept = 0, linetype = "dotted", color = "grey40") +
  
  theme_minimal(base_size = 15) +
  theme(
    panel.grid = element_blank(),
    axis.line = element_line(color = "black"),
    axis.text.y = element_text(color = "black", size = 14),  # Bigger Y-axis numbers
    axis.text.x = element_text(color = "black", size = 14),  # Bigger X-axis category names
    axis.title = element_text(color = "black", size = 16),   # Bigger axis titles
    legend.position = "none"
  ) +
  
  labs(x = "Pf infection status", y = "Total IgM+ Troph-iRBCs (% cells)")

p

ggsave(here("./figures/pf_infection_IgM.png"))
## Saving 7 x 5 in image

2 Final models

2.1 IgG

2.1.1 IgG+ iRBC

Proportion of IgG positive versus IgG negative infected red blood cells.

A binomial GLM is not applicable because of violated model assumptions.

model <- glm(
  cbind(raw_igGpos_iRBC,
        raw_igGneg_iRBC) ~ HBB.genotype,
  family = binomial, data = df)

summary(model)
## 
## Call:
## glm(formula = cbind(raw_igGpos_iRBC, raw_igGneg_iRBC) ~ HBB.genotype, 
##     family = binomial, data = df)
## 
## Coefficients:
##                   Estimate Std. Error  z value Pr(>|z|)    
## (Intercept)      -1.545394   0.009922 -155.762  < 2e-16 ***
## HBB.genotypeHbAC -0.020901   0.014188   -1.473    0.141    
## HBB.genotypeHbAS -0.219481   0.014104  -15.562  < 2e-16 ***
## HBB.genotypeHbCC  0.207452   0.026210    7.915 2.47e-15 ***
## HBB.genotypeHbSC -0.012781   0.021102   -0.606    0.545    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 8242.3  on 50  degrees of freedom
## Residual deviance: 7807.8  on 46  degrees of freedom
## AIC: 8237.9
## 
## Number of Fisher Scoring iterations: 4
# Simulate residuals
sim_res <- simulateResiduals(model, n = 1000)
plot(sim_res)

# Test for overdispersion
testDispersion(sim_res)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 8.0713, p-value < 2.2e-16
## alternative hypothesis: two.sided
# Test for zero-inflation
testZeroInflation(sim_res)

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = NaN, p-value = 1
## alternative hypothesis: two.sided
# plot residuals per group
plotResiduals(sim_res, df$HBB.genotype)

The beta-binomial GLM was used as the final model. The data was modelled using generalized linear model using the beta-binomial family and a single dispersion parameter to account for overdispersion.

model <- glmmTMB(
    cbind(raw_igGpos_iRBC,
          raw_igGneg_iRBC) ~ HBB.genotype,
  family = betabinomial(link = "logit"),
  data = df
)
summary(model)
##  Family: betabinomial  ( logit )
## Formula:          cbind(raw_igGpos_iRBC, raw_igGneg_iRBC) ~ HBB.genotype
## Data: df
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##     746.6     758.2    -367.3     734.6        45 
## 
## 
## Dispersion parameter for betabinomial family (): 27.1 
## 
## Conditional model:
##                  Estimate Std. Error z value Pr(>|z|)    
## (Intercept)      -1.52094    0.12918 -11.774   <2e-16 ***
## HBB.genotypeHbAC -0.10308    0.18413  -0.560    0.576    
## HBB.genotypeHbAS -0.24702    0.17918  -1.379    0.168    
## HBB.genotypeHbCC  0.24564    0.34265   0.717    0.473    
## HBB.genotypeHbSC -0.03087    0.27385  -0.113    0.910    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Simulate residuals
sim_res <- simulateResiduals(model, n = 1000)
plot(sim_res)

# Test for overdispersion
testDispersion(sim_res)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.81732, p-value = 0.42
## alternative hypothesis: two.sided
# Test for zero-inflation
testZeroInflation(sim_res)

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = NaN, p-value = 1
## alternative hypothesis: two.sided
# plot residuals per group
plotResiduals(sim_res, df$HBB.genotype)

The estimated odds ratios (using Wald Z-tests) for the different genotype contrasts are as follows:

# Compute estimated marginal means for the HBB.genotype factor
emm <- emmeans(model, ~ HBB.genotype)

# Perform all pairwise comparisons
pairs(emm, adjust="tukey", ratios=T, type="response", reverse=T, infer = c(TRUE, TRUE))
##  contrast    odds.ratio    SE  df asymp.LCL asymp.UCL null z.ratio p.value
##  HbAC / HbAA      0.902 0.166 Inf     0.546      1.49    1  -0.560  0.9807
##  HbAS / HbAA      0.781 0.140 Inf     0.479      1.27    1  -1.379  0.6414
##  HbAS / HbAC      0.866 0.158 Inf     0.527      1.42    1  -0.791  0.9333
##  HbCC / HbAA      1.278 0.438 Inf     0.502      3.26    1   0.717  0.9527
##  HbCC / HbAC      1.417 0.488 Inf     0.554      3.62    1   1.013  0.8494
##  HbCC / HbAS      1.637 0.559 Inf     0.645      4.16    1   1.442  0.6001
##  HbSC / HbAA      0.970 0.266 Inf     0.459      2.05    1  -0.113  1.0000
##  HbSC / HbAC      1.075 0.296 Inf     0.507      2.28    1   0.262  0.9990
##  HbSC / HbAS      1.241 0.338 Inf     0.590      2.61    1   0.793  0.9326
##  HbSC / HbCC      0.758 0.303 Inf     0.255      2.25    1  -0.692  0.9582
## 
## Confidence level used: 0.95 
## Conf-level adjustment: tukey method for comparing a family of 5 estimates 
## Intervals are back-transformed from the log odds ratio scale 
## P value adjustment: tukey method for comparing a family of 5 estimates 
## Tests are performed on the log odds ratio scale
# perform only comparisons against AA
AA = c(1,0,0,0,0)
AC = c(0,1,0,0,0)
AS = c(0,0,1,0,0)
CC = c(0,0,0,1,0)
SC = c(0,0,0,0,1)

contrast(emm, method = list("HbAC / HbAA" = AC - AA,
                            "HbAS / HbAA" = AS - AA,
                            "HbCC / HbAA" = CC - AA,
                            "HbSC / HbAA" = SC - AA),
         adjust = "bonferroni",
         type="response",
         infer = c(TRUE, TRUE)
         )
##  contrast    odds.ratio    SE  df asymp.LCL asymp.UCL null z.ratio p.value
##  HbAC / HbAA      0.902 0.166 Inf     0.570      1.43    1  -0.560  1.0000
##  HbAS / HbAA      0.781 0.140 Inf     0.499      1.22    1  -1.379  0.6720
##  HbCC / HbAA      1.278 0.438 Inf     0.543      3.01    1   0.717  1.0000
##  HbSC / HbAA      0.970 0.266 Inf     0.489      1.92    1  -0.113  1.0000
## 
## Confidence level used: 0.95 
## Conf-level adjustment: bonferroni method for 4 estimates 
## Intervals are back-transformed from the log odds ratio scale 
## P value adjustment: bonferroni method for 4 tests 
## Tests are performed on the log odds ratio scale

The above table not only shows the estimated odds ratio and its corresponding p-value, but also the confidence interval of this odds ratio (= asymp.LCL and asymp.UCL).

These p-values are adjusted for multiple testing using the Bonferroni method.

The estimates are already transformed from the log-scale, so they can be directly interpreted as odds ratios (OR).

The estimated SCR per group is, adjusted for multiple testing using the Bonferroni method:

confint(emm, type="response", adjust="bonferroni")

2.1.1.1 Bayesian estimates

brm_model <- brm(
  bf(
    raw_igGpos_iRBC | trials(raw_igGpos_iRBC + raw_igGneg_iRBC) ~ 
      HBB.genotype
  ), 
  family = beta_binomial(link = "logit"),
  data = df,
  iter = 5000,
  chains = 4,
  # control = list(adapt_delta = 0.95),
  seed = 42,
  cores = 8
)
## Compiling Stan program...
## Start sampling
summary(brm_model)
##  Family: beta_binomial 
##   Links: mu = logit 
## Formula: raw_igGpos_iRBC | trials(raw_igGpos_iRBC + raw_igGneg_iRBC) ~ HBB.genotype 
##    Data: df (Number of observations: 51) 
##   Draws: 4 chains, each with iter = 5000; warmup = 2500; thin = 1;
##          total post-warmup draws = 10000
## 
## Regression Coefficients:
##                  Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept           -1.52      0.14    -1.80    -1.25 1.00     9384     7216
## HBB.genotypeHbAC    -0.10      0.20    -0.49     0.29 1.00     9731     7487
## HBB.genotypeHbAS    -0.24      0.19    -0.62     0.14 1.00     9421     8149
## HBB.genotypeHbCC     0.18      0.40    -0.67     0.89 1.00    10567     6076
## HBB.genotypeHbSC    -0.06      0.31    -0.69     0.50 1.00    11056     6669
## 
## Further Distributional Parameters:
##     Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## phi    24.64      4.97    15.80    35.40 1.00     9980     7170
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
# Compute estimated marginal means for the HBB.genotype factor
emm <- emmeans(brm_model, ~ HBB.genotype)

# Perform all pairwise comparisons
pairs(emm, adjust="tukey", ratios=T, type="response", reverse=T, infer = c(TRUE, TRUE))
##  contrast    odds.ratio lower.HPD upper.HPD
##  HbAC / HbAA      0.903     0.586      1.30
##  HbAS / HbAA      0.783     0.499      1.10
##  HbAS / HbAC      0.868     0.571      1.24
##  HbCC / HbAA      1.230     0.417      2.27
##  HbCC / HbAC      1.368     0.469      2.55
##  HbCC / HbAS      1.578     0.563      2.94
##  HbSC / HbAA      0.955     0.466      1.60
##  HbSC / HbAC      1.057     0.509      1.78
##  HbSC / HbAS      1.221     0.564      2.02
##  HbSC / HbCC      0.776     0.238      1.74
## 
## Point estimate displayed: median 
## Results are back-transformed from the log odds ratio scale 
## HPD interval probability: 0.95
# perform only comparisons against AA
AA = c(1,0,0,0,0)
AC = c(0,1,0,0,0)
AS = c(0,0,1,0,0)
CC = c(0,0,0,1,0)
SC = c(0,0,0,0,1)

contrast(emm, method = list("HbAC / HbAA" = AC - AA,
                            "HbAS / HbAA" = AS - AA,
                            "HbCC / HbAA" = CC - AA,
                            "HbSC / HbAA" = SC - AA),
         adjust = "bonferroni",
         type="response",
         infer = c(TRUE, TRUE)
         )
##  contrast    odds.ratio lower.HPD upper.HPD
##  HbAC / HbAA      0.903     0.586      1.30
##  HbAS / HbAA      0.783     0.499      1.10
##  HbCC / HbAA      1.230     0.417      2.27
##  HbSC / HbAA      0.955     0.466      1.60
## 
## Point estimate displayed: median 
## Results are back-transformed from the log odds ratio scale 
## HPD interval probability: 0.95
confint(emm, type="response", adjust="bonferroni")

2.1.2 IgG+ asexual iRBC

The binomial model was not applicable.

model <- glm(
  cbind(raw_igGpos_asexual,
        raw_igGneg_asexual) ~ HBB.genotype,
  family = binomial, data = df)

summary(model)
## 
## Call:
## glm(formula = cbind(raw_igGpos_asexual, raw_igGneg_asexual) ~ 
##     HBB.genotype, family = binomial, data = df)
## 
## Coefficients:
##                  Estimate Std. Error  z value Pr(>|z|)    
## (Intercept)      -1.29328    0.01178 -109.789  < 2e-16 ***
## HBB.genotypeHbAC  0.01222    0.01680    0.727    0.467    
## HBB.genotypeHbAS -0.27538    0.01688  -16.310  < 2e-16 ***
## HBB.genotypeHbCC  0.23987    0.03097    7.745 9.54e-15 ***
## HBB.genotypeHbSC -0.10229    0.02549   -4.013 6.00e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 7977.4  on 50  degrees of freedom
## Residual deviance: 7482.8  on 46  degrees of freedom
## AIC: 7892.4
## 
## Number of Fisher Scoring iterations: 4
# Simulate residuals
sim_res <- simulateResiduals(model, n = 1000)
plot(sim_res)

# Test for overdispersion
testDispersion(sim_res)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 7.2991, p-value < 2.2e-16
## alternative hypothesis: two.sided
# Test for zero-inflation
testZeroInflation(sim_res)

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = NaN, p-value = 1
## alternative hypothesis: two.sided
# plot residuals per group
plotResiduals(sim_res, df$HBB.genotype)

The beta-binomial GLM was used as the final model. The data was modelled using generalized linear model using the beta-binomial family and a single dispersion parameter to account for overdispersion.

model <- glmmTMB(
    cbind(raw_igGpos_asexual,
          raw_igGneg_asexual) ~ HBB.genotype,
  family = betabinomial(link = "logit"),
  data = df
)
summary(model)
##  Family: betabinomial  ( logit )
## Formula:          cbind(raw_igGpos_asexual, raw_igGneg_asexual) ~ HBB.genotype
## Data: df
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##     725.7     737.3    -356.8     713.7        45 
## 
## 
## Dispersion parameter for betabinomial family (): 16.4 
## 
## Conditional model:
##                  Estimate Std. Error z value Pr(>|z|)    
## (Intercept)      -1.27711    0.15202  -8.401   <2e-16 ***
## HBB.genotypeHbAC -0.09105    0.21565  -0.422    0.673    
## HBB.genotypeHbAS -0.30215    0.21096  -1.432    0.152    
## HBB.genotypeHbCC  0.29470    0.40383   0.730    0.466    
## HBB.genotypeHbSC -0.09908    0.32566  -0.304    0.761    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Simulate residuals
sim_res <- simulateResiduals(model, n = 1000)
plot(sim_res)

# Test for overdispersion
testDispersion(sim_res)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.83947, p-value = 0.48
## alternative hypothesis: two.sided
# Test for zero-inflation
testZeroInflation(sim_res)

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = NaN, p-value = 1
## alternative hypothesis: two.sided
# plot residuals per group
plotResiduals(sim_res, df$HBB.genotype)

The estimated odds ratios (using Wald Z-tests) for the different genotype contrasts are as follows:

# Compute estimated marginal means for the HBB.genotype factor
emm <- emmeans(model, ~ HBB.genotype)

# Perform all pairwise comparisons
pairs(emm, adjust="tukey", ratios=T, type="response", reverse=T, infer = c(TRUE, TRUE))
##  contrast    odds.ratio    SE  df asymp.LCL asymp.UCL null z.ratio p.value
##  HbAC / HbAA      0.913 0.197 Inf     0.507      1.64    1  -0.422  0.9934
##  HbAS / HbAA      0.739 0.156 Inf     0.416      1.31    1  -1.432  0.6067
##  HbAS / HbAC      0.810 0.173 Inf     0.452      1.45    1  -0.989  0.8603
##  HbCC / HbAA      1.343 0.542 Inf     0.446      4.04    1   0.730  0.9496
##  HbCC / HbAC      1.471 0.596 Inf     0.487      4.44    1   0.952  0.8763
##  HbCC / HbAS      1.816 0.732 Inf     0.605      5.45    1   1.482  0.5743
##  HbSC / HbAA      0.906 0.295 Inf     0.373      2.20    1  -0.304  0.9981
##  HbSC / HbAC      0.992 0.325 Inf     0.406      2.42    1  -0.025  1.0000
##  HbSC / HbAS      1.225 0.397 Inf     0.506      2.97    1   0.627  0.9709
##  HbSC / HbCC      0.675 0.319 Inf     0.186      2.45    1  -0.833  0.9206
## 
## Confidence level used: 0.95 
## Conf-level adjustment: tukey method for comparing a family of 5 estimates 
## Intervals are back-transformed from the log odds ratio scale 
## P value adjustment: tukey method for comparing a family of 5 estimates 
## Tests are performed on the log odds ratio scale
# perform only comparisons against AA
AA = c(1,0,0,0,0)
AC = c(0,1,0,0,0)
AS = c(0,0,1,0,0)
CC = c(0,0,0,1,0)
SC = c(0,0,0,0,1)

contrast(emm, method = list("HbAC / HbAA" = AC - AA,
                            "HbAS / HbAA" = AS - AA,
                            "HbCC / HbAA" = CC - AA,
                            "HbSC / HbAA" = SC - AA),
         adjust = "bonferroni",
         type="response",
         infer = c(TRUE, TRUE)
         )
##  contrast    odds.ratio    SE  df asymp.LCL asymp.UCL null z.ratio p.value
##  HbAC / HbAA      0.913 0.197 Inf     0.533      1.56    1  -0.422  1.0000
##  HbAS / HbAA      0.739 0.156 Inf     0.436      1.25    1  -1.432  0.6083
##  HbCC / HbAA      1.343 0.542 Inf     0.490      3.68    1   0.730  1.0000
##  HbSC / HbAA      0.906 0.295 Inf     0.402      2.04    1  -0.304  1.0000
## 
## Confidence level used: 0.95 
## Conf-level adjustment: bonferroni method for 4 estimates 
## Intervals are back-transformed from the log odds ratio scale 
## P value adjustment: bonferroni method for 4 tests 
## Tests are performed on the log odds ratio scale

The above table not only shows the estimated odds ratio and its corresponding p-value, but also the confidence interval of this odds ratio (= asymp.LCL and asymp.UCL).

These p-values are adjusted for multiple testing using the Bonferroni method.

The estimates are already transformed from the log-scale, so they can be directly interpreted as odds ratios (OR).

The estimated SCR per group is, adjusted for multiple testing using the Bonferroni method:

confint(emm, type="response", adjust="bonferroni")

2.1.2.1 Bayesian estimates

brm_model <- brm(
  bf(
    raw_igGpos_asexual | trials(raw_igGpos_asexual + raw_igGneg_asexual) ~ 
      HBB.genotype
  ), 
  family = beta_binomial(link = "logit"),
  data = df,
  iter = 5000,
  chains = 4,
  # control = list(adapt_delta = 0.95),
  seed = 42,
  cores = 8
)
## Compiling Stan program...
## Start sampling
summary(brm_model)
##  Family: beta_binomial 
##   Links: mu = logit 
## Formula: raw_igGpos_asexual | trials(raw_igGpos_asexual + raw_igGneg_asexual) ~ HBB.genotype 
##    Data: df (Number of observations: 51) 
##   Draws: 4 chains, each with iter = 5000; warmup = 2500; thin = 1;
##          total post-warmup draws = 10000
## 
## Regression Coefficients:
##                  Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept           -1.27      0.16    -1.60    -0.96 1.00     8837     6807
## HBB.genotypeHbAC    -0.09      0.23    -0.53     0.36 1.00     8943     7844
## HBB.genotypeHbAS    -0.30      0.22    -0.74     0.13 1.00     8743     7089
## HBB.genotypeHbCC     0.23      0.48    -0.76     1.09 1.00     9564     6469
## HBB.genotypeHbSC    -0.14      0.35    -0.87     0.54 1.00     9960     7026
## 
## Further Distributional Parameters:
##     Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## phi    15.07      3.05     9.63    21.51 1.00    10702     7405
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
# Compute estimated marginal means for the HBB.genotype factor
emm <- emmeans(brm_model, ~ HBB.genotype)

# Perform all pairwise comparisons
pairs(emm, adjust="tukey", ratios=T, type="response", reverse=T, infer = c(TRUE, TRUE))
##  contrast    odds.ratio lower.HPD upper.HPD
##  HbAC / HbAA      0.915     0.545      1.37
##  HbAS / HbAA      0.741     0.455      1.10
##  HbAS / HbAC      0.807     0.473      1.19
##  HbCC / HbAA      1.296     0.340      2.71
##  HbCC / HbAC      1.416     0.385      2.96
##  HbCC / HbAS      1.749     0.451      3.61
##  HbSC / HbAA      0.888     0.370      1.59
##  HbSC / HbAC      0.967     0.413      1.76
##  HbSC / HbAS      1.192     0.511      2.14
##  HbSC / HbCC      0.683     0.151      1.76
## 
## Point estimate displayed: median 
## Results are back-transformed from the log odds ratio scale 
## HPD interval probability: 0.95
# perform only comparisons against AA
AA = c(1,0,0,0,0)
AC = c(0,1,0,0,0)
AS = c(0,0,1,0,0)
CC = c(0,0,0,1,0)
SC = c(0,0,0,0,1)

contrast(emm, method = list("HbAC / HbAA" = AC - AA,
                            "HbAS / HbAA" = AS - AA,
                            "HbCC / HbAA" = CC - AA,
                            "HbSC / HbAA" = SC - AA),
         adjust = "bonferroni",
         type="response",
         infer = c(TRUE, TRUE)
         )
##  contrast    odds.ratio lower.HPD upper.HPD
##  HbAC / HbAA      0.915     0.545      1.37
##  HbAS / HbAA      0.741     0.455      1.10
##  HbCC / HbAA      1.296     0.340      2.71
##  HbSC / HbAA      0.888     0.370      1.59
## 
## Point estimate displayed: median 
## Results are back-transformed from the log odds ratio scale 
## HPD interval probability: 0.95
confint(emm, type="response", adjust="bonferroni")

2.1.3 IgG+ sexual iRBC

The binomial model was not applicable.

model <- glm(
  cbind(raw_igGpos_sexual,
        raw_igGneg_sexual) ~ HBB.genotype,
  family = binomial, data = df)

summary(model)
## 
## Call:
## glm(formula = cbind(raw_igGpos_sexual, raw_igGneg_sexual) ~ HBB.genotype, 
##     family = binomial, data = df)
## 
## Coefficients:
##                  Estimate Std. Error z value Pr(>|z|)    
## (Intercept)      -3.72119    0.03966 -93.816  < 2e-16 ***
## HBB.genotypeHbAC -0.03870    0.05682  -0.681  0.49589    
## HBB.genotypeHbAS -0.15682    0.05613  -2.794  0.00521 ** 
## HBB.genotypeHbCC  0.22857    0.10188   2.243  0.02487 *  
## HBB.genotypeHbSC  0.02220    0.08398   0.264  0.79154    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 444.12  on 50  degrees of freedom
## Residual deviance: 426.00  on 46  degrees of freedom
## AIC: 715.45
## 
## Number of Fisher Scoring iterations: 4
# Simulate residuals
sim_res <- simulateResiduals(model, n = 1000)
plot(sim_res)

# Test for overdispersion
testDispersion(sim_res)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 4.6992, p-value < 2.2e-16
## alternative hypothesis: two.sided
# Test for zero-inflation
testZeroInflation(sim_res)

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = NaN, p-value = 1
## alternative hypothesis: two.sided
# plot residuals per group
plotResiduals(sim_res, df$HBB.genotype)

The beta-binomial GLM was used as the final model. The data was modelled using generalized linear model using the beta-binomial family and a single dispersion parameter to account for overdispersion.

model <- glmmTMB(
    cbind(raw_igGpos_sexual,
          raw_igGneg_sexual) ~ HBB.genotype,
  family = betabinomial(link = "logit"),
  data = df
)
summary(model)
##  Family: betabinomial  ( logit )
## Formula:          cbind(raw_igGpos_sexual, raw_igGneg_sexual) ~ HBB.genotype
## Data: df
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##     454.8     466.4    -221.4     442.8        45 
## 
## 
## Dispersion parameter for betabinomial family ():  235 
## 
## Conditional model:
##                  Estimate Std. Error z value Pr(>|z|)    
## (Intercept)      -3.69892    0.11681  -31.67   <2e-16 ***
## HBB.genotypeHbAC -0.11348    0.16838   -0.67    0.500    
## HBB.genotypeHbAS -0.15115    0.16111   -0.94    0.348    
## HBB.genotypeHbCC  0.27159    0.29402    0.92    0.356    
## HBB.genotypeHbSC  0.03024    0.24244    0.12    0.901    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Simulate residuals
sim_res <- simulateResiduals(model, n = 1000)
plot(sim_res)

# Test for overdispersion
testDispersion(sim_res)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.78507, p-value = 0.416
## alternative hypothesis: two.sided
# Test for zero-inflation
testZeroInflation(sim_res)

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = NaN, p-value = 1
## alternative hypothesis: two.sided
# plot residuals per group
plotResiduals(sim_res, df$HBB.genotype)

The estimated odds ratios (using Wald Z-tests) for the different genotype contrasts are as follows:

# Compute estimated marginal means for the HBB.genotype factor
emm <- emmeans(model, ~ HBB.genotype)

# Perform all pairwise comparisons
pairs(emm, adjust="tukey", ratios=T, type="response", reverse=T, infer = c(TRUE, TRUE))
##  contrast    odds.ratio    SE  df asymp.LCL asymp.UCL null z.ratio p.value
##  HbAC / HbAA      0.893 0.150 Inf     0.564      1.41    1  -0.674  0.9621
##  HbAS / HbAA      0.860 0.139 Inf     0.554      1.33    1  -0.938  0.8820
##  HbAS / HbAC      0.963 0.160 Inf     0.612      1.52    1  -0.227  0.9994
##  HbCC / HbAA      1.312 0.386 Inf     0.588      2.93    1   0.924  0.8878
##  HbCC / HbAC      1.470 0.436 Inf     0.654      3.30    1   1.297  0.6929
##  HbCC / HbAS      1.526 0.447 Inf     0.687      3.39    1   1.444  0.5992
##  HbSC / HbAA      1.031 0.250 Inf     0.532      2.00    1   0.125  0.9999
##  HbSC / HbAC      1.155 0.284 Inf     0.590      2.26    1   0.585  0.9774
##  HbSC / HbAS      1.199 0.289 Inf     0.621      2.31    1   0.753  0.9438
##  HbSC / HbCC      0.786 0.270 Inf     0.307      2.01    1  -0.701  0.9563
## 
## Confidence level used: 0.95 
## Conf-level adjustment: tukey method for comparing a family of 5 estimates 
## Intervals are back-transformed from the log odds ratio scale 
## P value adjustment: tukey method for comparing a family of 5 estimates 
## Tests are performed on the log odds ratio scale
# perform only comparisons against AA
AA = c(1,0,0,0,0)
AC = c(0,1,0,0,0)
AS = c(0,0,1,0,0)
CC = c(0,0,0,1,0)
SC = c(0,0,0,0,1)

contrast(emm, method = list("HbAC / HbAA" = AC - AA,
                            "HbAS / HbAA" = AS - AA,
                            "HbCC / HbAA" = CC - AA,
                            "HbSC / HbAA" = SC - AA),
         adjust = "bonferroni",
         type="response",
         infer = c(TRUE, TRUE)
         )
##  contrast    odds.ratio    SE  df asymp.LCL asymp.UCL null z.ratio p.value
##  HbAC / HbAA      0.893 0.150 Inf     0.586      1.36    1  -0.674  1.0000
##  HbAS / HbAA      0.860 0.139 Inf     0.575      1.29    1  -0.938  1.0000
##  HbCC / HbAA      1.312 0.386 Inf     0.630      2.73    1   0.924  1.0000
##  HbSC / HbAA      1.031 0.250 Inf     0.563      1.89    1   0.125  1.0000
## 
## Confidence level used: 0.95 
## Conf-level adjustment: bonferroni method for 4 estimates 
## Intervals are back-transformed from the log odds ratio scale 
## P value adjustment: bonferroni method for 4 tests 
## Tests are performed on the log odds ratio scale

The above table not only shows the estimated odds ratio and its corresponding p-value, but also the confidence interval of this odds ratio (= asymp.LCL and asymp.UCL).

These p-values are adjusted for multiple testing using the Bonferroni method.

The estimates are already transformed from the log-scale, so they can be directly interpreted as odds ratios (OR).

The estimated SCR per group is, adjusted for multiple testing using the Bonferroni method:

confint(emm, type="response", adjust="bonferroni")

2.1.3.1 Bayesian estimates

brm_model <- brm(
  bf(
    raw_igGpos_sexual | trials(raw_igGpos_sexual + raw_igGneg_sexual) ~ 
      HBB.genotype
  ), 
  family = beta_binomial(link = "logit"),
  data = df,
  iter = 5000,
  chains = 4,
  # control = list(adapt_delta = 0.95),
  seed = 42,
  cores = 8
)
## Compiling Stan program...
## Start sampling
summary(brm_model)
##  Family: beta_binomial 
##   Links: mu = logit 
## Formula: raw_igGpos_sexual | trials(raw_igGpos_sexual + raw_igGneg_sexual) ~ HBB.genotype 
##    Data: df (Number of observations: 51) 
##   Draws: 4 chains, each with iter = 5000; warmup = 2500; thin = 1;
##          total post-warmup draws = 10000
## 
## Regression Coefficients:
##                  Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept           -3.69      0.13    -3.95    -3.44 1.00     8685     7215
## HBB.genotypeHbAC    -0.12      0.19    -0.49     0.25 1.00     9199     7839
## HBB.genotypeHbAS    -0.15      0.18    -0.49     0.20 1.00     9227     7465
## HBB.genotypeHbCC     0.20      0.35    -0.59     0.81 1.00     9808     6006
## HBB.genotypeHbSC    -0.00      0.27    -0.59     0.50 1.00    10082     7296
## 
## Further Distributional Parameters:
##     Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## phi   192.71     45.35   115.79   291.89 1.00    11297     7268
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
# Compute estimated marginal means for the HBB.genotype factor
emm <- emmeans(brm_model, ~ HBB.genotype)

# Perform all pairwise comparisons
pairs(emm, adjust="tukey", ratios=T, type="response", reverse=T, infer = c(TRUE, TRUE))
##  contrast    odds.ratio lower.HPD upper.HPD
##  HbAC / HbAA      0.889     0.585      1.24
##  HbAS / HbAA      0.862     0.597      1.20
##  HbAS / HbAC      0.968     0.659      1.38
##  HbCC / HbAA      1.249     0.486      2.14
##  HbCC / HbAC      1.403     0.516      2.41
##  HbCC / HbAS      1.450     0.580      2.49
##  HbSC / HbAA      1.010     0.517      1.59
##  HbSC / HbAC      1.137     0.585      1.83
##  HbSC / HbAS      1.173     0.617      1.84
##  HbSC / HbCC      0.811     0.292      1.68
## 
## Point estimate displayed: median 
## Results are back-transformed from the log odds ratio scale 
## HPD interval probability: 0.95
# perform only comparisons against AA
AA = c(1,0,0,0,0)
AC = c(0,1,0,0,0)
AS = c(0,0,1,0,0)
CC = c(0,0,0,1,0)
SC = c(0,0,0,0,1)

contrast(emm, method = list("HbAC / HbAA" = AC - AA,
                            "HbAS / HbAA" = AS - AA,
                            "HbCC / HbAA" = CC - AA,
                            "HbSC / HbAA" = SC - AA),
         adjust = "bonferroni",
         type="response",
         infer = c(TRUE, TRUE)
         )
##  contrast    odds.ratio lower.HPD upper.HPD
##  HbAC / HbAA      0.889     0.585      1.24
##  HbAS / HbAA      0.862     0.597      1.20
##  HbCC / HbAA      1.249     0.486      2.14
##  HbSC / HbAA      1.010     0.517      1.59
## 
## Point estimate displayed: median 
## Results are back-transformed from the log odds ratio scale 
## HPD interval probability: 0.95
confint(emm, type="response", adjust="bonferroni")

2.2 IgM

2.2.1 IgM+ iRBC

Proportion of IgM positive versus IgM negative infected red blood cells.

A binomial GLM is not applicable because of violated model assumptions.

model <- glm(
  cbind(raw_igMpos_iRBC,
        raw_igMneg_iRBC) ~ HBB.genotype,
  family = binomial, data = df)

summary(model)
## 
## Call:
## glm(formula = cbind(raw_igMpos_iRBC, raw_igMneg_iRBC) ~ HBB.genotype, 
##     family = binomial, data = df)
## 
## Coefficients:
##                   Estimate Std. Error  z value Pr(>|z|)    
## (Intercept)      -3.355457   0.026214 -128.002  < 2e-16 ***
## HBB.genotypeHbAC  0.078928   0.036758    2.147 0.031773 *  
## HBB.genotypeHbAS  0.004638   0.035976    0.129 0.897430    
## HBB.genotypeHbCC -0.270520   0.081797   -3.307 0.000942 ***
## HBB.genotypeHbSC -0.610926   0.070615   -8.651  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 2713.6  on 50  degrees of freedom
## Residual deviance: 2587.5  on 46  degrees of freedom
## AIC: 2911.4
## 
## Number of Fisher Scoring iterations: 5
# Simulate residuals
sim_res <- simulateResiduals(model, n = 1000)
plot(sim_res)

# Test for overdispersion
testDispersion(sim_res)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 19.074, p-value < 2.2e-16
## alternative hypothesis: two.sided
# Test for zero-inflation
testZeroInflation(sim_res)

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = NaN, p-value = 1
## alternative hypothesis: two.sided
# plot residuals per group
plotResiduals(sim_res, df$HBB.genotype)

The beta-binomial GLM was used as the final model. The data was modelled using generalized linear model using the beta-binomial family and a single dispersion parameter to account for overdispersion.

model <- glmmTMB(
    cbind(raw_igMpos_iRBC,
          raw_igMneg_iRBC) ~ HBB.genotype,
  family = betabinomial(link = "logit"),
  data = df
)
summary(model)
##  Family: betabinomial  ( logit )
## Formula:          cbind(raw_igMpos_iRBC, raw_igMneg_iRBC) ~ HBB.genotype
## Data: df
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##     569.0     580.6    -278.5     557.0        45 
## 
## 
## Dispersion parameter for betabinomial family (): 66.2 
## 
## Conditional model:
##                  Estimate Std. Error z value Pr(>|z|)    
## (Intercept)      -3.32915    0.16831 -19.780   <2e-16 ***
## HBB.genotypeHbAC -0.06698    0.23416  -0.286    0.775    
## HBB.genotypeHbAS -0.03856    0.22235  -0.173    0.862    
## HBB.genotypeHbCC -0.06662    0.47091  -0.141    0.888    
## HBB.genotypeHbSC -0.43597    0.39062  -1.116    0.264    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Simulate residuals
sim_res <- simulateResiduals(model, n = 1000)
plot(sim_res)

# Test for overdispersion
testDispersion(sim_res)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 1.4807, p-value = 0.14
## alternative hypothesis: two.sided
# Test for zero-inflation
testZeroInflation(sim_res)

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 0, p-value = 1
## alternative hypothesis: two.sided
# plot residuals per group
plotResiduals(sim_res, df$HBB.genotype)

The estimated odds ratios (using Wald Z-tests) for the different genotype contrasts are as follows:

# Compute estimated marginal means for the HBB.genotype factor
emm <- emmeans(model, ~ HBB.genotype)

# Perform all pairwise comparisons
pairs(emm, adjust="tukey", ratios=T, type="response", reverse=T, infer = c(TRUE, TRUE))
##  contrast    odds.ratio    SE  df asymp.LCL asymp.UCL null z.ratio p.value
##  HbAC / HbAA      0.935 0.219 Inf     0.494      1.77    1  -0.286  0.9985
##  HbAS / HbAA      0.962 0.214 Inf     0.525      1.76    1  -0.173  0.9998
##  HbAS / HbAC      1.029 0.232 Inf     0.556      1.90    1   0.126  0.9999
##  HbCC / HbAA      0.936 0.441 Inf     0.259      3.38    1  -0.141  0.9999
##  HbCC / HbAC      1.000 0.472 Inf     0.276      3.63    1   0.001  1.0000
##  HbCC / HbAS      0.972 0.454 Inf     0.272      3.47    1  -0.060  1.0000
##  HbSC / HbAA      0.647 0.253 Inf     0.223      1.88    1  -1.116  0.7981
##  HbSC / HbAC      0.691 0.271 Inf     0.237      2.02    1  -0.941  0.8810
##  HbSC / HbAS      0.672 0.259 Inf     0.235      1.92    1  -1.031  0.8409
##  HbSC / HbCC      0.691 0.391 Inf     0.147      3.24    1  -0.652  0.9663
## 
## Confidence level used: 0.95 
## Conf-level adjustment: tukey method for comparing a family of 5 estimates 
## Intervals are back-transformed from the log odds ratio scale 
## P value adjustment: tukey method for comparing a family of 5 estimates 
## Tests are performed on the log odds ratio scale
# perform only comparisons against AA
AA = c(1,0,0,0,0)
AC = c(0,1,0,0,0)
AS = c(0,0,1,0,0)
CC = c(0,0,0,1,0)
SC = c(0,0,0,0,1)

contrast(emm, method = list("HbAC / HbAA" = AC - AA,
                            "HbAS / HbAA" = AS - AA,
                            "HbCC / HbAA" = CC - AA,
                            "HbSC / HbAA" = SC - AA),
         adjust = "bonferroni",
         type="response",
         infer = c(TRUE, TRUE)
         )
##  contrast    odds.ratio    SE  df asymp.LCL asymp.UCL null z.ratio p.value
##  HbAC / HbAA      0.935 0.219 Inf     0.521      1.68    1  -0.286  1.0000
##  HbAS / HbAA      0.962 0.214 Inf     0.552      1.68    1  -0.173  1.0000
##  HbCC / HbAA      0.936 0.441 Inf     0.289      3.03    1  -0.141  1.0000
##  HbSC / HbAA      0.647 0.253 Inf     0.244      1.72    1  -1.116  1.0000
## 
## Confidence level used: 0.95 
## Conf-level adjustment: bonferroni method for 4 estimates 
## Intervals are back-transformed from the log odds ratio scale 
## P value adjustment: bonferroni method for 4 tests 
## Tests are performed on the log odds ratio scale

The above table not only shows the estimated odds ratio and its corresponding p-value, but also the confidence interval of this odds ratio (= asymp.LCL and asymp.UCL).

These p-values are adjusted for multiple testing using the Bonferroni method.

The estimates are already transformed from the log-scale, so they can be directly interpreted as odds ratios (OR).

The estimated SCR per group is, adjusted for multiple testing using the Bonferroni method:

confint(emm, type="response", adjust="bonferroni")

2.2.1.1 Bayesian estimates

brm_model <- brm(
  bf(
    raw_igMpos_iRBC | trials(raw_igMpos_iRBC + raw_igMneg_iRBC) ~ 
      HBB.genotype
  ), 
  family = beta_binomial(link = "logit"),
  data = df,
  iter = 5000,
  chains = 4,
  # control = list(adapt_delta = 0.95),
  seed = 42,
  cores = 8
)
## Compiling Stan program...
## Start sampling
summary(brm_model)
##  Family: beta_binomial 
##   Links: mu = logit 
## Formula: raw_igMpos_iRBC | trials(raw_igMpos_iRBC + raw_igMneg_iRBC) ~ HBB.genotype 
##    Data: df (Number of observations: 51) 
##   Draws: 4 chains, each with iter = 5000; warmup = 2500; thin = 1;
##          total post-warmup draws = 10000
## 
## Regression Coefficients:
##                  Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept           -3.32      0.18    -3.68    -2.98 1.00     8818     7200
## HBB.genotypeHbAC    -0.06      0.25    -0.54     0.43 1.00     9089     7851
## HBB.genotypeHbAS    -0.03      0.23    -0.48     0.44 1.00     9147     7077
## HBB.genotypeHbCC    -0.23      0.57    -1.53     0.72 1.00     9201     5612
## HBB.genotypeHbSC    -0.51      0.44    -1.47     0.28 1.00     9252     6355
## 
## Further Distributional Parameters:
##     Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## phi    60.07     12.87    37.31    87.21 1.00     9920     6906
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
# Compute estimated marginal means for the HBB.genotype factor
emm <- emmeans(brm_model, ~ HBB.genotype)

# Perform all pairwise comparisons
pairs(emm, adjust="tukey", ratios=T, type="response", reverse=T, infer = c(TRUE, TRUE))
##  contrast    odds.ratio lower.HPD upper.HPD
##  HbAC / HbAA      0.939    0.5227      1.44
##  HbAS / HbAA      0.971    0.5785      1.48
##  HbAS / HbAC      1.031    0.5993      1.55
##  HbCC / HbAA      0.851    0.1254      1.85
##  HbCC / HbAC      0.895    0.1284      1.97
##  HbCC / HbAS      0.868    0.1642      1.92
##  HbSC / HbAA      0.622    0.1656      1.21
##  HbSC / HbAC      0.662    0.1874      1.29
##  HbSC / HbAS      0.640    0.1845      1.23
##  HbSC / HbCC      0.737    0.0844      2.45
## 
## Point estimate displayed: median 
## Results are back-transformed from the log odds ratio scale 
## HPD interval probability: 0.95
# perform only comparisons against AA
AA = c(1,0,0,0,0)
AC = c(0,1,0,0,0)
AS = c(0,0,1,0,0)
CC = c(0,0,0,1,0)
SC = c(0,0,0,0,1)

contrast(emm, method = list("HbAC / HbAA" = AC - AA,
                            "HbAS / HbAA" = AS - AA,
                            "HbCC / HbAA" = CC - AA,
                            "HbSC / HbAA" = SC - AA),
         adjust = "bonferroni",
         type="response",
         infer = c(TRUE, TRUE)
         )
##  contrast    odds.ratio lower.HPD upper.HPD
##  HbAC / HbAA      0.939     0.523      1.44
##  HbAS / HbAA      0.971     0.578      1.48
##  HbCC / HbAA      0.851     0.125      1.85
##  HbSC / HbAA      0.622     0.166      1.21
## 
## Point estimate displayed: median 
## Results are back-transformed from the log odds ratio scale 
## HPD interval probability: 0.95
confint(emm, type="response", adjust="bonferroni")

2.2.2 IgM+ asexual iRBC

The binomial model was not applicable.

model <- glm(
  cbind(raw_igMpos_asexual,
        raw_igMneg_asexual) ~ HBB.genotype,
  family = binomial, data = df)

summary(model)
## 
## Call:
## glm(formula = cbind(raw_igMpos_asexual, raw_igMneg_asexual) ~ 
##     HBB.genotype, family = binomial, data = df)
## 
## Coefficients:
##                   Estimate Std. Error  z value Pr(>|z|)    
## (Intercept)      -3.372290   0.026616 -126.700  < 2e-16 ***
## HBB.genotypeHbAC  0.068161   0.037417    1.822  0.06851 .  
## HBB.genotypeHbAS -0.001644   0.036614   -0.045  0.96419    
## HBB.genotypeHbCC -0.272719   0.083005   -3.286  0.00102 ** 
## HBB.genotypeHbSC -0.590718   0.071156   -8.302  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 2613.4  on 50  degrees of freedom
## Residual deviance: 2500.1  on 46  degrees of freedom
## AIC: 2822.2
## 
## Number of Fisher Scoring iterations: 5
# Simulate residuals
sim_res <- simulateResiduals(model, n = 1000)
plot(sim_res)

# Test for overdispersion
testDispersion(sim_res)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 19.528, p-value < 2.2e-16
## alternative hypothesis: two.sided
# Test for zero-inflation
testZeroInflation(sim_res)

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = NaN, p-value = 1
## alternative hypothesis: two.sided
# plot residuals per group
plotResiduals(sim_res, df$HBB.genotype)

The beta-binomial GLM was used as the final model. The data was modelled using generalized linear model using the beta-binomial family and a single dispersion parameter to account for overdispersion.

model <- glmmTMB(
    cbind(raw_igMpos_asexual,
          raw_igMneg_asexual) ~ HBB.genotype,
  family = betabinomial(link = "logit"),
  data = df
)
summary(model)
##  Family: betabinomial  ( logit )
## Formula:          cbind(raw_igMpos_asexual, raw_igMneg_asexual) ~ HBB.genotype
## Data: df
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##     565.0     576.6    -276.5     553.0        45 
## 
## 
## Dispersion parameter for betabinomial family (): 68.2 
## 
## Conditional model:
##                  Estimate Std. Error z value Pr(>|z|)    
## (Intercept)      -3.34908    0.16756 -19.988   <2e-16 ***
## HBB.genotypeHbAC -0.07570    0.23357  -0.324    0.746    
## HBB.genotypeHbAS -0.03852    0.22143  -0.174    0.862    
## HBB.genotypeHbCC -0.06949    0.46938  -0.148    0.882    
## HBB.genotypeHbSC -0.42438    0.38807  -1.094    0.274    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Simulate residuals
sim_res <- simulateResiduals(model, n = 1000)
plot(sim_res)

# Test for overdispersion
testDispersion(sim_res)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 1.533, p-value = 0.08
## alternative hypothesis: two.sided
# Test for zero-inflation
testZeroInflation(sim_res)

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 0, p-value = 1
## alternative hypothesis: two.sided
# plot residuals per group
plotResiduals(sim_res, df$HBB.genotype)

The estimated odds ratios (using Wald Z-tests) for the different genotype contrasts are as follows:

# Compute estimated marginal means for the HBB.genotype factor
emm <- emmeans(model, ~ HBB.genotype)

# Perform all pairwise comparisons
pairs(emm, adjust="tukey", ratios=T, type="response", reverse=T, infer = c(TRUE, TRUE))
##  contrast    odds.ratio    SE  df asymp.LCL asymp.UCL null z.ratio p.value
##  HbAC / HbAA      0.927 0.217 Inf     0.490      1.75    1  -0.324  0.9976
##  HbAS / HbAA      0.962 0.213 Inf     0.526      1.76    1  -0.174  0.9998
##  HbAS / HbAC      1.038 0.233 Inf     0.562      1.92    1   0.165  0.9998
##  HbCC / HbAA      0.933 0.438 Inf     0.259      3.36    1  -0.148  0.9999
##  HbCC / HbAC      1.006 0.474 Inf     0.278      3.64    1   0.013  1.0000
##  HbCC / HbAS      0.969 0.451 Inf     0.273      3.45    1  -0.067  1.0000
##  HbSC / HbAA      0.654 0.254 Inf     0.227      1.89    1  -1.094  0.8100
##  HbSC / HbAC      0.706 0.275 Inf     0.244      2.04    1  -0.894  0.8991
##  HbSC / HbAS      0.680 0.260 Inf     0.239      1.93    1  -1.008  0.8518
##  HbSC / HbCC      0.701 0.395 Inf     0.151      3.26    1  -0.630  0.9703
## 
## Confidence level used: 0.95 
## Conf-level adjustment: tukey method for comparing a family of 5 estimates 
## Intervals are back-transformed from the log odds ratio scale 
## P value adjustment: tukey method for comparing a family of 5 estimates 
## Tests are performed on the log odds ratio scale
# perform only comparisons against AA
AA = c(1,0,0,0,0)
AC = c(0,1,0,0,0)
AS = c(0,0,1,0,0)
CC = c(0,0,0,1,0)
SC = c(0,0,0,0,1)

contrast(emm, method = list("HbAC / HbAA" = AC - AA,
                            "HbAS / HbAA" = AS - AA,
                            "HbCC / HbAA" = CC - AA,
                            "HbSC / HbAA" = SC - AA),
         adjust = "bonferroni",
         type="response",
         infer = c(TRUE, TRUE)
         )
##  contrast    odds.ratio    SE  df asymp.LCL asymp.UCL null z.ratio p.value
##  HbAC / HbAA      0.927 0.217 Inf     0.517      1.66    1  -0.324  1.0000
##  HbAS / HbAA      0.962 0.213 Inf     0.553      1.67    1  -0.174  1.0000
##  HbCC / HbAA      0.933 0.438 Inf     0.289      3.01    1  -0.148  1.0000
##  HbSC / HbAA      0.654 0.254 Inf     0.248      1.72    1  -1.094  1.0000
## 
## Confidence level used: 0.95 
## Conf-level adjustment: bonferroni method for 4 estimates 
## Intervals are back-transformed from the log odds ratio scale 
## P value adjustment: bonferroni method for 4 tests 
## Tests are performed on the log odds ratio scale

The above table not only shows the estimated odds ratio and its corresponding p-value, but also the confidence interval of this odds ratio (= asymp.LCL and asymp.UCL).

These p-values are adjusted for multiple testing using the Bonferroni method.

The estimates are already transformed from the log-scale, so they can be directly interpreted as odds ratios (OR).

The estimated SCR per group is, adjusted for multiple testing using the Bonferroni method:

confint(emm, type="response", adjust="bonferroni")

2.2.2.1 Bayesian estimates

brm_model <- brm(
  bf(
    raw_igMpos_asexual | trials(raw_igMpos_asexual + raw_igMneg_asexual) ~ 
      HBB.genotype
  ), 
  family = beta_binomial(link = "logit"),
  data = df,
  iter = 5000,
  chains = 4,
  # control = list(adapt_delta = 0.95),
  seed = 42,
  cores = 8
)
## Compiling Stan program...
## Start sampling
summary(brm_model)
##  Family: beta_binomial 
##   Links: mu = logit 
## Formula: raw_igMpos_asexual | trials(raw_igMpos_asexual + raw_igMneg_asexual) ~ HBB.genotype 
##    Data: df (Number of observations: 51) 
##   Draws: 4 chains, each with iter = 5000; warmup = 2500; thin = 1;
##          total post-warmup draws = 10000
## 
## Regression Coefficients:
##                  Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept           -3.34      0.18    -3.71    -3.00 1.00     8254     6656
## HBB.genotypeHbAC    -0.07      0.25    -0.56     0.42 1.00     8688     7179
## HBB.genotypeHbAS    -0.03      0.23    -0.47     0.44 1.00     8979     7311
## HBB.genotypeHbCC    -0.24      0.58    -1.60     0.70 1.00     9310     5117
## HBB.genotypeHbSC    -0.49      0.43    -1.44     0.28 1.00     9819     5964
## 
## Further Distributional Parameters:
##     Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## phi    61.74     13.42    38.01    91.19 1.00     9051     6501
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
# Compute estimated marginal means for the HBB.genotype factor
emm <- emmeans(brm_model, ~ HBB.genotype)

# Perform all pairwise comparisons
pairs(emm, adjust="tukey", ratios=T, type="response", reverse=T, infer = c(TRUE, TRUE))
##  contrast    odds.ratio lower.HPD upper.HPD
##  HbAC / HbAA      0.932    0.5367      1.45
##  HbAS / HbAA      0.968    0.5748      1.46
##  HbAS / HbAC      1.043    0.6174      1.57
##  HbCC / HbAA      0.842    0.1035      1.79
##  HbCC / HbAC      0.899    0.1303      1.97
##  HbCC / HbAS      0.866    0.1012      1.84
##  HbSC / HbAA      0.635    0.1707      1.20
##  HbSC / HbAC      0.680    0.2016      1.30
##  HbSC / HbAS      0.655    0.1688      1.20
##  HbSC / HbCC      0.759    0.0971      2.45
## 
## Point estimate displayed: median 
## Results are back-transformed from the log odds ratio scale 
## HPD interval probability: 0.95
# perform only comparisons against AA
AA = c(1,0,0,0,0)
AC = c(0,1,0,0,0)
AS = c(0,0,1,0,0)
CC = c(0,0,0,1,0)
SC = c(0,0,0,0,1)

contrast(emm, method = list("HbAC / HbAA" = AC - AA,
                            "HbAS / HbAA" = AS - AA,
                            "HbCC / HbAA" = CC - AA,
                            "HbSC / HbAA" = SC - AA),
         adjust = "bonferroni",
         type="response",
         infer = c(TRUE, TRUE)
         )
##  contrast    odds.ratio lower.HPD upper.HPD
##  HbAC / HbAA      0.932     0.537      1.45
##  HbAS / HbAA      0.968     0.575      1.46
##  HbCC / HbAA      0.842     0.103      1.79
##  HbSC / HbAA      0.635     0.171      1.20
## 
## Point estimate displayed: median 
## Results are back-transformed from the log odds ratio scale 
## HPD interval probability: 0.95
confint(emm, type="response", adjust="bonferroni")

2.2.3 IgM+ sexual iRBC

The binomial model was not applicable.

model <- glm(
  cbind(raw_igMpos_sexual,
        raw_igMneg_sexual) ~ HBB.genotype,
  family = binomial, data = df)

summary(model)
## 
## Call:
## glm(formula = cbind(raw_igMpos_sexual, raw_igMneg_sexual) ~ HBB.genotype, 
##     family = binomial, data = df)
## 
## Coefficients:
##                   Estimate Std. Error z value Pr(>|z|)    
## (Intercept)        -7.1787     0.2183 -32.884   <2e-16 ***
## HBB.genotypeHbAC    0.4958     0.2777   1.785   0.0742 .  
## HBB.genotypeHbAS    0.1753     0.2867   0.611   0.5409    
## HBB.genotypeHbCC   -1.1068     1.0237  -1.081   0.2796    
## HBB.genotypeHbSC  -17.0827  1278.5349  -0.013   0.9893    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 166.65  on 50  degrees of freedom
## Residual deviance: 146.32  on 46  degrees of freedom
## AIC: 227.54
## 
## Number of Fisher Scoring iterations: 16
# Simulate residuals
sim_res <- simulateResiduals(model, n = 1000)
plot(sim_res)

# Test for overdispersion
testDispersion(sim_res)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 2.9299, p-value < 2.2e-16
## alternative hypothesis: two.sided
# Test for zero-inflation
testZeroInflation(sim_res)

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 2.0402, p-value < 2.2e-16
## alternative hypothesis: two.sided
# plot residuals per group
plotResiduals(sim_res, df$HBB.genotype)

The beta-binomial GLM was used as the final model. The data was modelled using generalized linear model using the beta-binomial family and a single dispersion parameter to account for overdispersion.

model <- glmmTMB(
    cbind(raw_igMpos_sexual,
          raw_igMneg_sexual) ~ HBB.genotype,
  family = betabinomial(link = "logit"),
  data = df
)
summary(model)
##  Family: betabinomial  ( logit )
## Formula:          cbind(raw_igMpos_sexual, raw_igMneg_sexual) ~ HBB.genotype
## Data: df
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##     180.0     191.5     -84.0     168.0        45 
## 
## 
## Dispersion parameter for betabinomial family ():  573 
## 
## Conditional model:
##                    Estimate Std. Error z value Pr(>|z|)    
## (Intercept)      -7.007e+00  3.688e-01 -18.999   <2e-16 ***
## HBB.genotypeHbAC  1.406e-01  4.868e-01   0.289    0.773    
## HBB.genotypeHbAS  2.055e-02  4.631e-01   0.044    0.965    
## HBB.genotypeHbCC -4.401e-01  1.057e+00  -0.416    0.677    
## HBB.genotypeHbSC -1.974e+01  1.107e+04  -0.002    0.999    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Simulate residuals
sim_res <- simulateResiduals(model, n = 1000)
plot(sim_res)

# Test for overdispersion
testDispersion(sim_res)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.80213, p-value = 0.83
## alternative hypothesis: two.sided
# Test for zero-inflation
testZeroInflation(sim_res)

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 1.0213, p-value = 0.92
## alternative hypothesis: two.sided
# plot residuals per group
plotResiduals(sim_res, df$HBB.genotype)

The estimated odds ratios (using Wald Z-tests) for the different genotype contrasts are as follows:

# Compute estimated marginal means for the HBB.genotype factor
emm <- emmeans(model, ~ HBB.genotype)

# Perform all pairwise comparisons
pairs(emm, adjust="tukey", ratios=T, type="response", reverse=T, infer = c(TRUE, TRUE))
##  contrast    odds.ratio       SE  df asymp.LCL asymp.UCL null z.ratio p.value
##  HbAC / HbAA      1.151 5.60e-01 Inf    0.3051         4    1   0.289  0.9985
##  HbAS / HbAA      1.021 4.73e-01 Inf    0.2886         4    1   0.044  1.0000
##  HbAS / HbAC      0.887 4.17e-01 Inf    0.2462         3    1  -0.256  0.9991
##  HbCC / HbAA      0.644 6.81e-01 Inf    0.0360        12    1  -0.416  0.9937
##  HbCC / HbAC      0.560 5.95e-01 Inf    0.0307        10    1  -0.546  0.9825
##  HbCC / HbAS      0.631 6.63e-01 Inf    0.0358        11    1  -0.438  0.9924
##  HbSC / HbAA      0.000 2.95e-05 Inf    0.0000       Inf    1  -0.002  1.0000
##  HbSC / HbAC      0.000 2.56e-05 Inf    0.0000       Inf    1  -0.002  1.0000
##  HbSC / HbAS      0.000 2.89e-05 Inf    0.0000       Inf    1  -0.002  1.0000
##  HbSC / HbCC      0.000 4.57e-05 Inf    0.0000       Inf    1  -0.002  1.0000
## 
## Confidence level used: 0.95 
## Conf-level adjustment: tukey method for comparing a family of 5 estimates 
## Intervals are back-transformed from the log odds ratio scale 
## P value adjustment: tukey method for comparing a family of 5 estimates 
## Tests are performed on the log odds ratio scale
# perform only comparisons against AA
AA = c(1,0,0,0,0)
AC = c(0,1,0,0,0)
AS = c(0,0,1,0,0)
CC = c(0,0,0,1,0)
SC = c(0,0,0,0,1)

contrast(emm, method = list("HbAC / HbAA" = AC - AA,
                            "HbAS / HbAA" = AS - AA,
                            "HbCC / HbAA" = CC - AA,
                            "HbSC / HbAA" = SC - AA),
         adjust = "bonferroni",
         type="response",
         infer = c(TRUE, TRUE)
         )
##  contrast    odds.ratio       SE  df asymp.LCL asymp.UCL null z.ratio p.value
##  HbAC / HbAA      1.151 5.60e-01 Inf    0.3412         4    1   0.289  1.0000
##  HbAS / HbAA      1.021 4.73e-01 Inf    0.3210         3    1   0.044  1.0000
##  HbCC / HbAA      0.644 6.81e-01 Inf    0.0459         9    1  -0.416  1.0000
##  HbSC / HbAA      0.000 2.95e-05 Inf    0.0000       Inf    1  -0.002  1.0000
## 
## Confidence level used: 0.95 
## Conf-level adjustment: bonferroni method for 4 estimates 
## Intervals are back-transformed from the log odds ratio scale 
## P value adjustment: bonferroni method for 4 tests 
## Tests are performed on the log odds ratio scale

The above table not only shows the estimated odds ratio and its corresponding p-value, but also the confidence interval of this odds ratio (= asymp.LCL and asymp.UCL).

These p-values are adjusted for multiple testing using the Bonferroni method.

The estimates are already transformed from the log-scale, so they can be directly interpreted as odds ratios (OR).

The estimated SCR per group is, adjusted for multiple testing using the Bonferroni method:

confint(emm, type="response", adjust="bonferroni")

2.2.3.1 Bayesian estimates

brm_model <- brm(
  bf(
    raw_igMpos_sexual | trials(raw_igMpos_sexual + raw_igMneg_sexual) ~ 
      HBB.genotype
  ), 
  family = beta_binomial(link = "logit"),
  data = df,
  iter = 5000,
  chains = 4,
  # control = list(adapt_delta = 0.95),
  seed = 42,
  cores = 8
)
## Compiling Stan program...
## Start sampling
## Warning: There were 37 divergent transitions after warmup. See
## https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## to find out why this is a problem and how to eliminate them.
## Warning: Examine the pairs() plot to diagnose sampling problems
summary(brm_model)
## Warning: There were 37 divergent transitions after warmup. Increasing
## adapt_delta above 0.8 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
##  Family: beta_binomial 
##   Links: mu = logit 
## Formula: raw_igMpos_sexual | trials(raw_igMpos_sexual + raw_igMneg_sexual) ~ HBB.genotype 
##    Data: df (Number of observations: 51) 
##   Draws: 4 chains, each with iter = 5000; warmup = 2500; thin = 1;
##          total post-warmup draws = 10000
## 
## Regression Coefficients:
##                  Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept           -6.72      0.41    -7.55    -5.94 1.00     3695     4338
## HBB.genotypeHbAC     0.04      0.51    -0.98     1.04 1.00     3673     3306
## HBB.genotypeHbAS    -0.01      0.49    -0.97     0.98 1.00     3907     3915
## HBB.genotypeHbCC    -0.85      1.29    -4.01     1.16 1.00     3809     3003
## HBB.genotypeHbSC   -57.62     72.78  -270.09    -2.42 1.00     1057      711
## 
## Further Distributional Parameters:
##     Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## phi   294.73    114.34   112.63   561.06 1.00     4388     4357
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
# Compute estimated marginal means for the HBB.genotype factor
emm <- emmeans(brm_model, ~ HBB.genotype)

# Perform all pairwise comparisons
pairs(emm, adjust="tukey", ratios=T, type="response", reverse=T, infer = c(TRUE, TRUE))
##  contrast    odds.ratio lower.HPD upper.HPD
##  HbAC / HbAA      1.043  0.277451    2.4730
##  HbAS / HbAA      0.984  0.262777    2.2540
##  HbAS / HbAC      0.950  0.239007    2.2416
##  HbCC / HbAA      0.520  0.000906    2.4482
##  HbCC / HbAC      0.495  0.000608    2.3643
##  HbCC / HbAS      0.520  0.000762    2.4449
##  HbSC / HbAA      0.000  0.000000    0.0272
##  HbSC / HbAC      0.000  0.000000    0.0264
##  HbSC / HbAS      0.000  0.000000    0.0277
##  HbSC / HbCC      0.000  0.000000    0.0702
## 
## Point estimate displayed: median 
## Results are back-transformed from the log odds ratio scale 
## HPD interval probability: 0.95
# perform only comparisons against AA
AA = c(1,0,0,0,0)
AC = c(0,1,0,0,0)
AS = c(0,0,1,0,0)
CC = c(0,0,0,1,0)
SC = c(0,0,0,0,1)

contrast(emm, method = list("HbAC / HbAA" = AC - AA,
                            "HbAS / HbAA" = AS - AA,
                            "HbCC / HbAA" = CC - AA,
                            "HbSC / HbAA" = SC - AA),
         adjust = "bonferroni",
         type="response",
         infer = c(TRUE, TRUE)
         )
##  contrast    odds.ratio lower.HPD upper.HPD
##  HbAC / HbAA      1.043  0.277451    2.4730
##  HbAS / HbAA      0.984  0.262777    2.2540
##  HbCC / HbAA      0.520  0.000906    2.4482
##  HbSC / HbAA      0.000  0.000000    0.0272
## 
## Point estimate displayed: median 
## Results are back-transformed from the log odds ratio scale 
## HPD interval probability: 0.95
confint(emm, type="response", adjust="bonferroni")

2.3 Pf infection

2.3.1 IgG+

Proportion of IgG positive versus IgG negative infected red blood cells for active vs non-active infections:

A binomial GLM is not applicable because of violated model assumptions.

model <- glm(
  cbind(raw_igGpos_iRBC,
        raw_igGneg_iRBC) ~ Pf_infection,
  family = binomial, data = df)

summary(model)
## 
## Call:
## glm(formula = cbind(raw_igGpos_iRBC, raw_igGneg_iRBC) ~ Pf_infection, 
##     family = binomial, data = df)
## 
## Coefficients:
##                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)     -1.96465    0.01127 -174.35   <2e-16 ***
## Pf_infectionYES  0.48053    0.01284   37.42   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 8242.3  on 50  degrees of freedom
## Residual deviance: 6752.8  on 49  degrees of freedom
## AIC: 7176.8
## 
## Number of Fisher Scoring iterations: 4
# Simulate residuals
sim_res <- simulateResiduals(model, n = 1000)
plot(sim_res)

# Test for overdispersion
testDispersion(sim_res)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 4.6219, p-value < 2.2e-16
## alternative hypothesis: two.sided
# Test for zero-inflation
testZeroInflation(sim_res)

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = NaN, p-value = 1
## alternative hypothesis: two.sided
# plot residuals per group
plotResiduals(sim_res, df$Pf_infection)

The beta-binomial GLM was used as the final model. The data was modelled using generalized linear model using the beta-binomial family and a single dispersion parameter to account for overdispersion.

model <- glmmTMB(
  cbind(raw_igGpos_iRBC,
        raw_igGneg_iRBC) ~ Pf_infection,
  family = betabinomial(link = "logit"),
  data = df
)
summary(model)
##  Family: betabinomial  ( logit )
## Formula:          cbind(raw_igGpos_iRBC, raw_igGneg_iRBC) ~ Pf_infection
## Data: df
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##     728.2     734.0    -361.1     722.2        48 
## 
## 
## Dispersion parameter for betabinomial family (): 34.8 
## 
## Conditional model:
##                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)      -2.0615     0.1327 -15.532  < 2e-16 ***
## Pf_infectionYES   0.5944     0.1490   3.989 6.63e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Simulate residuals
sim_res <- simulateResiduals(model, n = 1000)
plot(sim_res)

# Test for overdispersion
testDispersion(sim_res)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.87804, p-value = 0.568
## alternative hypothesis: two.sided
# Test for zero-inflation
testZeroInflation(sim_res)

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = NaN, p-value = 1
## alternative hypothesis: two.sided
# plot residuals per group
plotResiduals(sim_res, df$Pf_infection)

The estimated odds ratios (using Wald Z-tests) for the different genotype contrasts are as follows:

# Compute estimated marginal means for the Pf_infection factor
emm <- emmeans(model, ~ Pf_infection)

# Perform all pairwise comparisons
pairs(emm, adjust="bonferroni", ratios=T, type="response", reverse=T, infer = c(TRUE, TRUE))
##  contrast odds.ratio   SE  df asymp.LCL asymp.UCL null z.ratio p.value
##  YES / NO       1.81 0.27 Inf      1.35      2.43    1   3.989  0.0001
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log odds ratio scale 
## Tests are performed on the log odds ratio scale

The above table not only shows the estimated odds ratio and its corresponding p-value, but also the confidence interval of this odds ratio (= asymp.LCL and asymp.UCL).

These p-values are adjusted for multiple testing using the Bonferroni method.

The estimates are already transformed from the log-scale, so they can be directly interpreted as odds ratios (OR).

The estimated SCR per group is, adjusted for multiple testing using the Bonferroni method:

confint(emm, type="response", adjust="bonferroni")

2.3.1.1 Bayesian estimates

brm_model <- brm(
  bf(
    raw_igGpos_iRBC | trials(raw_igGpos_iRBC + raw_igGneg_iRBC) ~ 
      Pf_infection
  ), 
  family = beta_binomial(link = "logit"),
  data = df,
  iter = 5000,
  chains = 4,
  # control = list(adapt_delta = 0.95),
  seed = 42,
  cores = 8
)
## Compiling Stan program...
## Start sampling
summary(brm_model)
##  Family: beta_binomial 
##   Links: mu = logit 
## Formula: raw_igGpos_iRBC | trials(raw_igGpos_iRBC + raw_igGneg_iRBC) ~ Pf_infection 
##    Data: df (Number of observations: 51) 
##   Draws: 4 chains, each with iter = 5000; warmup = 2500; thin = 1;
##          total post-warmup draws = 10000
## 
## Regression Coefficients:
##                 Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept          -2.06      0.14    -2.34    -1.80 1.00     6802     6563
## Pf_infectionYES     0.60      0.16     0.30     0.91 1.00     7860     6694
## 
## Further Distributional Parameters:
##     Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## phi    33.18      6.71    21.18    47.73 1.00     7335     6405
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
# Compute estimated marginal means for the Pf_infection factor
emm <- emmeans(brm_model, ~ Pf_infection)

# Perform all pairwise comparisons
pairs(emm, adjust="bonferroni", ratios=T, type="response", reverse=T, infer = c(TRUE, TRUE))
##  contrast odds.ratio lower.HPD upper.HPD
##  YES / NO       1.82      1.31      2.43
## 
## Point estimate displayed: median 
## Results are back-transformed from the log odds ratio scale 
## HPD interval probability: 0.95
confint(emm, type="response", adjust="bonferroni")

2.3.1.2 Model comparison

model <- glmmTMB(
  cbind(raw_igGpos_iRBC,
        raw_igGneg_iRBC) ~ Pf_infection,
  family = betabinomial(link = "logit"),
  data = df
)
model_extended <- glmmTMB(
  cbind(raw_igGpos_iRBC,
        raw_igGneg_iRBC) ~ Pf_infection + Gender + Village + Age_group,
  family = betabinomial(link = "logit"),
  data = df
)
anova(model, model_extended)

A more complex model is not preferred. Moreover, the directionality and magnitude of the effect remains consistent.

summary(model_extended)
##  Family: betabinomial  ( logit )
## Formula:          
## cbind(raw_igGpos_iRBC, raw_igGneg_iRBC) ~ Pf_infection + Gender +  
##     Village + Age_group
## Data: df
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##     728.5     743.9    -356.2     712.5        43 
## 
## 
## Dispersion parameter for betabinomial family (): 42.2 
## 
## Conditional model:
##                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)     -1.90038    0.15024 -12.649  < 2e-16 ***
## Pf_infectionYES  0.66309    0.14878   4.457 8.32e-06 ***
## Gender1         -0.03191    0.11609  -0.275  0.78339    
## Village2        -0.33062    0.13737  -2.407  0.01609 *  
## Village3        -0.23119    0.20133  -1.148  0.25083    
## Village4        -0.52382    0.19835  -2.641  0.00827 ** 
## Age_group2      -0.03899    0.12991  -0.300  0.76409    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Simulate residuals
sim_res <- simulateResiduals(model_extended, n = 1000)
plot(sim_res)

# Test for overdispersion
testDispersion(sim_res)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.84876, p-value = 0.462
## alternative hypothesis: two.sided
# Test for zero-inflation
testZeroInflation(sim_res)

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = NaN, p-value = 1
## alternative hypothesis: two.sided
# plot residuals per group
plotResiduals(sim_res, df$Pf_infection)

2.3.2 IgM+

Proportion of IgM positive versus IgM negative infected red blood cells for active vs non-active infections:

A binomial GLM is not applicable because of violated model assumptions.

model <- glm(
  cbind(raw_igMpos_iRBC,
        raw_igMneg_iRBC) ~ Pf_infection,
  family = binomial, data = df)

summary(model)
## 
## Call:
## glm(formula = cbind(raw_igMpos_iRBC, raw_igMneg_iRBC) ~ Pf_infection, 
##     family = binomial, data = df)
## 
## Coefficients:
##                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)     -3.87580    0.03309 -117.11   <2e-16 ***
## Pf_infectionYES  0.64817    0.03660   17.71   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 2713.6  on 50  degrees of freedom
## Residual deviance: 2358.5  on 49  degrees of freedom
## AIC: 2676.5
## 
## Number of Fisher Scoring iterations: 5
# Simulate residuals
sim_res <- simulateResiduals(model, n = 1000)
plot(sim_res)

# Test for overdispersion
testDispersion(sim_res)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 8.3282, p-value < 2.2e-16
## alternative hypothesis: two.sided
# Test for zero-inflation
testZeroInflation(sim_res)

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = NaN, p-value = 1
## alternative hypothesis: two.sided
# plot residuals per group
plotResiduals(sim_res, df$Pf_infection)

The beta-binomial GLM was used as the final model. The data was modelled using generalized linear model using the beta-binomial family and a single dispersion parameter to account for overdispersion.

model <- glmmTMB(
  cbind(raw_igMpos_iRBC,
        raw_igMneg_iRBC) ~ Pf_infection,
  family = betabinomial(link = "logit"),
  data = df
)
summary(model)
##  Family: betabinomial  ( logit )
## Formula:          cbind(raw_igMpos_iRBC, raw_igMneg_iRBC) ~ Pf_infection
## Data: df
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##     560.1     565.9    -277.1     554.1        48 
## 
## 
## Dispersion parameter for betabinomial family (): 70.5 
## 
## Conditional model:
##                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)      -3.6959     0.1830 -20.192   <2e-16 ***
## Pf_infectionYES   0.4089     0.2016   2.029   0.0425 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Simulate residuals
sim_res <- simulateResiduals(model, n = 1000)
plot(sim_res)

# Test for overdispersion
testDispersion(sim_res)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 1.5207, p-value = 0.126
## alternative hypothesis: two.sided
# Test for zero-inflation
testZeroInflation(sim_res)

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 0, p-value = 1
## alternative hypothesis: two.sided
# plot residuals per group
plotResiduals(sim_res, df$Pf_infection)

The estimated odds ratios (using Wald Z-tests) for the different genotype contrasts are as follows:

# Compute estimated marginal means for the Pf_infection factor
emm <- emmeans(model, ~ Pf_infection)

# Perform all pairwise comparisons
pairs(emm, adjust="bonferroni", ratios=T, type="response", reverse=T, infer = c(TRUE, TRUE))
##  contrast odds.ratio    SE  df asymp.LCL asymp.UCL null z.ratio p.value
##  YES / NO       1.51 0.303 Inf      1.01      2.23    1   2.029  0.0425
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log odds ratio scale 
## Tests are performed on the log odds ratio scale

The above table not only shows the estimated odds ratio and its corresponding p-value, but also the confidence interval of this odds ratio (= asymp.LCL and asymp.UCL).

These p-values are adjusted for multiple testing using the Bonferroni method.

The estimates are already transformed from the log-scale, so they can be directly interpreted as odds ratios (OR).

The estimated SCR per group is, adjusted for multiple testing using the Bonferroni method:

confint(emm, type="response", adjust="bonferroni")

2.3.2.1 Bayesian estimates

brm_model <- brm(
  bf(
    raw_igMpos_iRBC | trials(raw_igMpos_iRBC + raw_igMneg_iRBC) ~ 
      Pf_infection
  ), 
  family = beta_binomial(link = "logit"),
  data = df,
  iter = 5000,
  chains = 4,
  # control = list(adapt_delta = 0.95),
  seed = 42,
  cores = 8
)
## Compiling Stan program...
## Start sampling
summary(brm_model)
##  Family: beta_binomial 
##   Links: mu = logit 
## Formula: raw_igMpos_iRBC | trials(raw_igMpos_iRBC + raw_igMneg_iRBC) ~ Pf_infection 
##    Data: df (Number of observations: 51) 
##   Draws: 4 chains, each with iter = 5000; warmup = 2500; thin = 1;
##          total post-warmup draws = 10000
## 
## Regression Coefficients:
##                 Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept          -3.69      0.19    -4.07    -3.34 1.00     6274     6111
## Pf_infectionYES     0.42      0.21     0.03     0.83 1.00     7019     6525
## 
## Further Distributional Parameters:
##     Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## phi    66.14     14.05    41.43    95.92 1.00     5906     5889
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
# Compute estimated marginal means for the Pf_infection factor
emm <- emmeans(brm_model, ~ Pf_infection)

# Perform all pairwise comparisons
pairs(emm, adjust="bonferroni", ratios=T, type="response", reverse=T, infer = c(TRUE, TRUE))
##  contrast odds.ratio lower.HPD upper.HPD
##  YES / NO       1.51     0.963       2.2
## 
## Point estimate displayed: median 
## Results are back-transformed from the log odds ratio scale 
## HPD interval probability: 0.95
confint(emm, type="response", adjust="bonferroni")

3 Software used

grateful::cite_packages(output = "table", out.dir = here(), cite.tidyverse = TRUE)